'''
Created on Aug 3, 2009

@author: mkiyer
'''

import sys
from math import log10

# Solexa->Sanger quality conversion table
conv_table = []
for x in xrange(-64, 64):
    conv_table.append(chr(int(33 + 10 * log10(1 + 10**(x/10.0))/log10(10)+0.499)))
    
#my @conv_table;
#for (-64..64) {
#  $conv_table[$_+64] = chr(int(33 + 10*log(1+10**($_/10.0))/log(10)+.499));
#}

def convert_scarf_to_fastq(fin, fout):
    '''
    convert an old solexa sequence file in the format: 
    SEQ_ID:SEQUENCE:PHRED_SCORES
    
    to standard fastq format:
    @SEQ_ID
    SEQUENCE
    +
    PHRED_SCORES

    Example:
    in: "USI-EAS6_5:4:1:210:465:AAATTCCTGGACACATACACTCTCCCAAGACTA:27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 26 27 27 27 27 27 27 27 21 18 26 21 27 27 12"
    out:
    @USI-EAS6_5:4:1:210:465
    AAATTCCTGGACACATACACTCTCCCAAGACTA
    +
    27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 27 26 27 27 27 27 27 27 27 21 18 26 21 27 27 12
    '''    
    for thisline in fin:    
        fields = thisline.strip().split(':')
        fout.write('@%s\n%s\n+\n' % ('_'.join(fields[0:5]), fields[5]))
        quals = map(int, fields[6].split())
        qualstr = ''.join([conv_table[x+64] for x in quals])
        fout.write('%s\n' % qualstr)

def convert_elandexport_to_fastq(fin, fout):
    for thisline in fin:        
        fields = thisline.strip().split('\t')        
        name = ':'.join(fields[0:6])
        seq = fields[8]
        qual = fields[9]
        fout.write('@%s\n%s\n+%s\n%s\n' % (name, seq, name, qual))
        
#sub scarf2std {
#  while (<>) {
#    my @t = split(':', $_);
#    my $name = join('_', @t[0..4]);
#    print "\@$name\n$t[5]\n+\n";
#    my $qual = '';
#    @t = split(/\s/, $t[6]);
#    $qual .= $conv_table[$_+64] for (@t);
#    print "$qual\n";
#  }
#}

if __name__ == '__main__':
    cmd = sys.argv[1]
    fin = open(sys.argv[2])
    fout = sys.stdout
    if cmd == 'scarf':
        convert_scarf_to_fastq(fin, fout)
    elif cmd == 'export':
        convert_elandexport_to_fastq(fin, fout)